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In this paper, we present a numerical solution to an ordinary differential 
equation of a fractional order in one-dimensional space. The solution to this 
equation can describe a steady state of the process of anomalous diffusion. 
The process arises from interactions within complex and non-homogeneous 
background. We present a numerical method which is based on the finite 
differences method. We consider a boundary value problem (Dirichlet condi- 
tions) for an equation with the Riesz-Feller fractional derivative. In the final 
part of this paper, same simulation results are shown. We present an example 
of non-linear temperature profiles in nanotubes which can be approximated 
by a solution to the fractional differential equation. 
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1 Introduction 

Anomalous diffusion is a phenomenon strongly connected with interac- 
tions within complex and non-homogeneous background. This phenomenon 
is observed in transport of a fluid in porous materials, in chaotic heat baths, 
amorphous semiconductors, particle dynamics inside a polymer network, two- 
dimensional rotating flow and also in econophysics. The phenomenon of 
anomalous diffusion deviates from the standard diffusion behaviour. In op- 
posite to the standard diffusion where a linear form in the mean square 
displacement (x 2 (t)) ~ k±t of the diffusing particle over time occurs, the 



anomalous diffusion is characterized by a non-linear one {x 2 (£)) ~ /c 7 t 7 , for 
7 G (0, 2]. In this phenomenon, there may exist a dependence {x 2 (t)) — > oo, 
which is characterized by occurrence of rare but extremely large jumps of 
diffusing particles - well-known as the Levy motion or the Levy flights. An 
ordinary diffusion follows Gaussian statistics and Fick's second law for find- 
ing the running process at time t whereas anomalous diffusion follows non- 
Gaussian statistic or can be interpreted as the Levy stable densities. 

Many authors proposed models which base on linear and non-linear forms 
of differential equations. Such models can simulate anomalous diffusion 
but they do not reflect its real behaviour. Several authors (Carpinteri and 
Mainardi, 1997; Hilfer, 2000; Gorenflo and Mainardi, 1998; Metzler and 
Klafter, 2000) apply fractional calculus to the modelling of this type of diffu- 
sion. This means that time and spatial derivatives in the classical diffusion 
equation are replaced by fractional ones. In comparison to derivatives of 
the integer order, which depend on the local behaviour of the function, the 
derivatives of the fractional order accumulate the whole history of this func- 
tion. 

In our previeus works (Ciesielski and Leszczynski, 2003, 2005), we pre- 
sented a solution to a partial differential equation of the fractional order with 
the time fractional operator and the space ordinary operator, respectively. 
Those solutions were based on the Finite Difference Method (FDM) and are 
called the Fractional FDM (FFDM). 



2 Mathematical background 

In this paper we consider an ordinary differential equation of fractional 
order in the following form 

d „T(x) = iGl (2.1) 



d \x iv 

where T(x) is a field variable (i.e. field temperature), (d a /d\x\g)T(x) is 
the Riesz- Feller fractional operator (Metzer and Klafter, 2000; Samko et al, 
1993), a is the real order of this operator and 9 is the skewness parame- 
ter. According to (Gorenflo and Mainardi, 1998), the Riesz- Feller fractional 
operator is defined as 

' T(x) = X D«T (x) = - [c L (a, 6) ^D^T (x) + c R (a, 9) x D a +(y0 T (x)] 



d \x 



(2.2) 
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for < « < 2, a 7^ 1 where 

(j \ m 
[-oo/r Q T(x)] (2.3) 



for m G N, m — 1 < a < m, and the coefficients ci(a,9), cr{ol,6) (for 
< a < 2, a ^ 1, |#| < min (a, 2 — a)) are defined as 

(a -6) Ti {a + 9)u 
sm sm 

c L (a,6) = . ( 2 , c R (a,9) = . , 2 . (2.4) 

sm {an) sm (Q7r) 

The fractional integral operators of the order a: _oo I"T (x) and X I^T (x) 
are defined as the left- and right-hand of Weyl's fractional integrals (Carpin- 
teri and Mainardi, 1997; Oldham and Spaner, 1974; Podlubny, 1999; Samko 
et al, 1993) whose definitions are 



T(a) J-oo (x-0 1 ' 



.^T (x) = — j - ^_J £ (2.5) 



x 00 w r(a) J x (e-x) 1 - 

Considering Eq. (2.1), we obtain the steady state of the classical diffusion 
equation for a = 2, i.e. the heat transfer equation. If a = 1, and the 
parameter of skewness 9 admits extreme values in (|2.4|) . the steady state of a 
transport equation is noted. Therefore we assume variations of the parameter 
a within the range < a < 2. Analysing behaviour of the parameter a < 2 
in Eq. (2.1), we found some combination between transport and propagation 
processes in steady states. 

In this work we will consider equation (|2.1j) in one dimensional domain 
Q : L < x < R with the boundary-value conditions of the first kind (Dirichlet 
conditions) as 

X = L k T T?m = 9L ( 2 - 6 ) 

x = R : T (R) = g R 

3 Numerical method 

According to the finite difference method (Hoffman, 1992; Majchrzak and 
Mochnacki, 1996), we consider a discrete from of equation ()2.1j) in space. The 
problem of solving equation ()2.1j) lies in a properly approximation of Riesz- 
Feller derivative (|2.2j) by a numerical scheme. 
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3.1 Approximation of the Riesz-Feller derivative 

In order to develop a discrete form of operator (|2.2j) . we introduce a 
homogenous spatial grid — oo < . . . < £i_ 2 < < %i < ^i+i < < 
. . . < oo with the step h = Xk — Xk-i- We denote a value of the function T 
in the point Xk as = T (x^), for fceZ. For numerical integration scheme, 
we assumed the trapezoidal rule and we used various weighted numerical 
differential schemes for the first and second derivatives, respectively. The 
method of determination of the discrete form of this operator was described 
in detail in work by Ciesielski (2005) and its final form is the following 



(a,9) 



(3.1) 



where coefficients w^' e \ 



for < a < 1, have the form 



w 



(a,0) 



2T (2 

((1*1 + 2) 



a) 



l-a 



Ai + 



+ 1) 1 - Q (2-3A 1 ) 
+ Ifcl 1 "" (3X 1 - 4) + (|A;| - I) 1 "" (2 - Ai)) c L 
(3 1 -°A 1 + 2 l ~ a (2 - 3Ai) + 3Ai - 4) c L + X iCr 
(2 1 - a X l - 3Ai + 2) (c L + c R ) 
(3 1 ~ a X 1 + 2 l - a (2 - 3Ai) + 3Ai - 4) c R + X iCl 



{{k + 2f- a A x + (*+!) 



l-Q 



+fc 1 ~ Q (3Ai 
and for 1 < a < 2, we obtain 



(2-3A0 
4 ) + (A;-l) 1 - a (2 -AO) c R 



(3.2) 



for k < -2 
for k = —1 
for fc = 
for = 1 

for fc > 2 



-1 



2r (3 - a) 



-1) 
-6) 



lf- a (4A 2 - 6) + 
2 " a (4A 2 - 2) + 



6)c L + 



+ 2) 2 ^ Q (2 - A 2 ) 
+ |fc| 2 - Q (6-6A 2 )- 
+ {\k\-2f- a {-X 2 ))c L 
(3 2 - a (2 - A 2 ) + 2 2 - a (4A 2 - 6) - 6A 2 

+ (2 - A 2 ) c R 
(2 2 -"(2-A 2 )+4A 2 -6) (c L + c fl ) 
(3 2 - a (2 - A 2 ) + 2 2 - a (4A 2 - 6) - 6A 2 + 6) c R 

+ (2 - A 2 ) c L 
((k + 2f- a (2 - A 2 ) + (k + lf- a (4A 2 - 6) + 



+k 2 ~ a (6 — 6A 2 ) + (k — l) 
+ (k- 2f- a (-A 2 )) 



2-a 



(4A 2 - 2) + 



for k < -2 

for k = —1 
for = 

for k = 1 



for fc > 2 



(3.3) 
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Assuming a = 2 and 9 = and cl (2, 0) = c# (2, 0) = —1/2, we obtain 

for k < -2 
for k = —1 

for fc = (3.4) 
for k = 1 
for jfe > 2 

These coeeficients are identical as in the well known central difference scheme 
for the second derivative. 

The authors did not find exact values of the approximating coefficients 
in literature. When a = 1, the Riesz- Feller operator is singular, hence the 
problem occurs. 

3.2 Fractional FDM 

Having the discretization of the Riesz-Feller derivative in space done, in 
this subsection we describe the finite difference method for equation (|2.1jl . 
Here we restrict the solution to one dimensional space in comparison with the 
standard diffusion equation where the discretization of the second derivative 
in space can approximate the central difference of yhe second order. The 
differences appear in the setting of boundary conditions. 

In the scheme of FDM, we replaced equation (J2.1)) by the following for- 
mula 

^ £ T ^' e) = ( 3 - 5 ) 

k=— oo 

Here, for the unbounded domain we are obligated to solve a system of alge- 
braic equations with an infinite dimension. 

3.2.1 Boundary value problem 

Present numerical scheme (13. 5|) with the included unbounded domain 
— oo < x < oo has no practical implementations in computer simulations. 

L R 



Figure 1: Spatial distribution of grid nodes 

Now we show solution to this problem on the finite domain Q : L < x < R 
with boundary conditions (j2.6|) . We divide domain Q into iV sub-domains 



w 



(2,0) 




1 
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with h = (R — L)/N. Figure 1 shows a modified spatial grid. Here we can 
observe additional 'virtual' points in the grid placed outside of the domain 
Q. In order to introduce the Dirichlet boundary conditions, we proposed a 
treatment which is based on the assumption that values of the function T in 
outside points are identical as values in the boundary nodes x or x^ 

T( Xk )^\l^\=^ f r *<» (3.6) 
[ T (xn) = gR for k > N 

On the basis of above considerations, we modify expression (j3.1j) for the 
discretization of the Riesz-Feller derivative. Thus we have 

Xi D%T ( Xi ) « ^ ( £ T ^' 9) + + 9rSr { n'-I\ (3-7) 

for i = 1, . . . , N — 1, where 

( a ,e) v-^ 1 ( a ,e) _ J c L (a, 6) -Vj for < a < 1 



s 4 = 2^ w k 



k=-oo (. L V"' ^ 3 



cl (a, 0) ■ rj for 1 < a < 2 



(3.8) 



0,0) (a,0) _ J c R (a, 0) • rj for < a < 1 



and 



fc=i+1 ^ c R (a, 0) • Tj for 1 < a < 2 



U + 2) 1 -" Ax + (j + l) 1 -" (2 - 2AQ + j 1 - (Ai - 2) 
2r (2 - a) 

(j + 2) 2 " a (2 - A 2 ) + (j + lf- a (3A 2 - 4) + j 2 ~ a (2 - 3A 2 ) 



(3.9) 



2r (3 - a) 



{ 3 -lf- a \ 



2 



2r (3 - a) 



Putting expression (|3.7J) to equation ()2.1|) . we obtain the following finite 
difference scheme 



N-i 



£ + 9LS L i + 9r s Rn—I = 0, for i = 0, . . . , iV 



k=—i 



T = g L (3.10) 
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The above scheme described by expression (|3.10|) can be written in a matrix 
form as 



where 



A 
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T = 


1 











a_i 


a Q 


fli 


«2 


a_ 2 


a_i 


a 
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a_4 


0-3 


a_2 


a_i 



B 



(3.11) 



B 
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-iV+5 • • 
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fli 


«2 
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-N+2 


a-JV+3 a - 


-7V+4 • • 


a_i 


fl 
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L, 61,62, 


63, 


64, • ■ 


■ , 6aT-2, 6jv- 
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OiV-3 
OjV-4 
a 7V-5 
a 7V-6 





OiV-3 
0^-3 





OiV-2 
0<N-2 



(3.12) 



(a,0) . (a,0) 



for j 



for j 



(3.13) 



iV-1 



and T = [T , 7\, T 2 , ■ ■ ■ , T N ] is the vector of unknown values of the function 
T. 

We can observe that boundary conditions influence to all values of the 
function in every node. In opposite to the second derivative over space, 
which is approximated locally, the characteristic feature of the Riesz-Feller 
and other fractional derivatives is the dependence on values of all domain 
points. For a = 2 and 9 = 0, our scheme is identical as with the well known 
and used central difference scheme in space (Hoffman, 1992; Majchrzak and 
Mochnacki, 1996). 

The skewness parameter 9 has significant influence on the solution. For 
a — > 1 + and 9 — > ±1 + one can obtain the classical ordinary differential 
equation. 



4 Simulation results 

In this section, we present results of calculation. In all presented sim- 
ulations we assumed < x < 1. Figure 2 shows temperature profiles 
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over space with boundary conditions Ql = 2, = 1 for different values 
a = {0.1,0.5.0.75,1.01,1.25,1.5,1.75,2} and = 0. Figure 3 presents an- 
other example of the solution in which we assumed a = 1.01 and different 
values of the skewness parameter = {0, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99}. 
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Figure 2: Spatial solution to equation (2.1) for different values of a and = 




Figure 3: Spatial solution to equation (2.1) for a steady value of the param- 
eter a = 1.01 and different values of 9 

The last example illustrates the heat transport in nanotubes. Figure 4 
presents experimental data prerformed by Zhang and Li (2005) and results of 



8 



our numerical solution of equation (|2.1j) . For such a comparison we assumed 
a = 0.35 and 9 = —0.055 to best fit the experimental data. It should be 
noted that the temperature profile inside the nanotube deviates from the 
profile obtained by solving the standard heat transfer equation. 
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Figure 4: Comparison of temperature profiles of nanotubes measured by 
Zhang and Li (2005) and obtained by numerical solution to equation (2.1) 



5 Conclusions 

Summing up, we proposed a fractional finite difference method for the 
fractional diffusion equation with the Riesz- Feller fractional derivative which 
is an extension of the standard diffusion. We analysed a linear case of the 
steady state of the anomalous diffusion equation, and in the future we will 
work on non-linear cases. We obtained the FDM scheme which generalises 
classical scheme of FDM for the diffusion equation. Moreover, for a = 2, 
our solution is equivalent with that obtained by the classical finite difference 
method. 

Analysing plots included in this work, we can see that in the case a = 2, 
solution of Eq. (2.1) is a linear function. In the other cases, when a < 2, 
solutions are non-linear functions. When we analyse the probability density 
function generated by fractional diffusion equation for a < 2, we observe 
a long tail of distribution. In this way we can simulate same rare and ex- 
treme events which are characterised by arbitrary very large values of particle 
jumps. 
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Analysing the changes in the skewness parameter 0, we observed interest- 
ing behaviour in the solution. For a — > 1 + and for 6 — > ±1 + , we obtained the 
steady state of the first order wave equation. For 6 e (0, 1) (with restrictions 
to the order a), we generated non-symmetric solutions. It should be noted 
that we can good approximate the temperature profile inside the nanotubes 
using solution of Eq. (2.1). 
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